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We consider a general method for computing the sum of positive Lya- 
punov exponents for moderately dense gases. This method is based upon hi- 
erarchy techniques used previously to derive the generalized Boltzmann equa- 
tion for the time dependent spatial and velocity distribution functions for such 
systems. We extend the variables in the generalized Boltzmann equation to 
include a new set of quantities that describe the separation of trajectories in 
phase space needed for a calculation of the Lyapunov exponents. The method 
described here is especially suitable for calculating the sum of all of the posi- 
tive Lyapunov exponents for the system, and may be applied to equilibrium 
as well as non-equilibrium situations. For low densities we obtain an extended 
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Boltzmann equation, from which, under a simplifying approximation, we re- 
cover the sum of positive Lyapunov exponents for hard disk and hard sphere 
systems, obtained before by a simpler method. In addition we indicate how to 
improve these results by avoiding the simplifying approximation. The restric- 
tion to hard sphere systems in (^-dimensions is made to keep the somewhat 
complicated formalism as clear as possible, but the method can be easily gen- 
eralized to apply to gases of particles that interact with strong short range 
forces. 

I. INTRODUCTION 

The ergodic and mixing properties of classical <i-dimensional hard sphere systems has 
been the subject of a great deal of work over the past several decades. Stimulated by the 
progress made by Sinai and co-workers toward proving the ergodicity of billiard type systems 
in the 1960's and onward JI|, others have extended these results in a number of directions, 
including a recent proof of the ergodic behavior of some systems of hard spheres in 2 and 3 
dimensions by Szasz and Simanyi f2|. Moreover, numerical studies of the chaotic behavior 
of hard sphere systems have recently been carried out by Dellago, Posch and Hoover who 
computed the Lyapunov spectrum, and the Kolmogorov-Sinai (KS) entropy, h^s m the 
case that the spheres are in equilibrium or in non-equilibrium steady states 0. For closed 
systems, according to Pesin's theorem, the KS entropy is equal to the sum of the positive 
Lyapunov exponents . 

For a simple model system, the Lorentz gas at low densities, it has been possible to 
calculate analytically the Lyapunov spectrum and to show that the results so obtained are 
in good agreement with computer results 0. This model considers the motion of a point 
particle moving in a random array of fixed scatterers, which are usually taken to be hard 
spheres in d dimensions. It was shown that standard kinetic theory methods for treating the 
Lorentz gas could be extended to include a systematic calculation of the Lyapunov spectrum. 
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Methods used included simple free path analyses as well as an extension of the Lorentz- 
Boltzmann transport equation to include new variables which describe the separation of 
trajectories in the phase space of the moving particle. Thus the question arises as to whether 
it is possible to extend the Boltzmann transport equation and its generalization to higher 
densities, as well as other related methods in kinetic theory, to compute properties of the 
Lyapunov spectrum of a gas where all of the particles move. Some results in this direction 
have recently been obtained. 

A theoretical treatment of the Lyapunov spectrum and hxs has recently been undertaken 
for hard disks or spheres at low density, na d << 1, where n is the number density of 
the particles and a is their diameter. Van Zon and Van Beijeren made an approximate 
calculation of the largest Lyapunov exponent for a dilute gas of hard disks and obtained 
results that are in very good agreement with the computer simulations of Dellago ||. In 
addition van Beijeren, Dorfman, Posch, and Dellago presented an approximate calculation 
of hxs/N, the Kolmogorov-Sinai entropy per particle, where iV >> 1 is the number of 
particles in the system, for dilute gases of particles with short range forces in equilibrium 
||. For hard sphere and hard disk systems their calculations led to expressions for the KS 
entropy per particle of the form 

h KS /N = au[-\nn + b + 0(n)}. (1) 

Here v is the average collision frequency per particle in the gas; n = na 2 for d = 2, and 
n = rnra 3 for d = 3 is the reduced density of scatterers; and a, b are numerical constants 
depending on the number of spatial dimensions of the system. The values of a are in 
excellent agreement with computer simulations but the theoretical values of b differed from 
the computer results by a factor of about 3 for d = 3 and about 6 for d = 2. 

The purpose of this paper is to provide a systematic method for computing hxs/N, or in 
more general circumstances, the sum of all the positive Lyapunov exponents, for a system of 
iV particles interacting with short range forces. In principle, this method should give exact 
values for the coefficients a and b above, as well as higher density corrections to these results. 
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It employs BBGKY hierarchy techniques for calculating the KS entropy per particle. Very 
similar hierarchy methods have been used in the kinetic theory of moderately dense gases 
to derive the Boltzmann transport equation for f(r,v,t), and to generalize this equation to 
moderately dense gases ||10|| . Here / is the one particle distribution function which gives the 
probability of finding a particle in the gas at the spatial point r, with velocity v at time 
t. Here we show that it is possible to include new variables in the distribution function 
which describe the separation of trajectory bundles in the phase space of N particles. We 
can then extend and generalize the Boltzmann equation so as to provide a method for 
computing Hks f° r dilute and moderately dense gases in equilibrium and non-equilibrium 
steady states. In order to simplify the discussion we will consider hard-sphere systems in d 
dimensions and briefly indicate at various points how these results can be extended to more 
general potentials. 

This paper is organized as follows: In Section II we describe the separation of two 
trajectories in phase space in terms of a set of d x d matrices whose elements have the 
dimension of a time, and which we call the radius of curvature (ROC) matrices. We then 
show that the sum of the positive Lyapunov exponents can, for a hard sphere system in d 
dimensions, be written in terms of the ROC matrices and a reduced two-particle distribution 
function. In Section III we describe the time dependence of these matrices due to free motion 
of the particles and due to binary collisions. After showing that the time dependence of the 
ROC matrices can be expressed in terms of free streaming and binary collision operators, we 
show that the reduced distribution functions defined in Section II satisfy a set of hierarchy 
equations, similar to those appearing in the usual BBGKY hierarchy equations, but modified 
to include the ROC matrices. In Section IV we consider the low density case, and show that 
the sum of positive Lyapunov exponents at low densities can be computed in terms of a single 
particle distribution function which satisfies an extended Boltzmann transport equation. In 
Section V we show how the previous results of Van Beijeren, Dorfman, Posch, and Dellago |9| 
for the KS entropy per particle, for an equilibrium system at low densities, can be obtained 
from a partial solution of this extended Boltzmann equation for equilibrium systems, and 
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argue that an important contribution is neglected in this approximate solution. Due to 
the complicated method that is needed to compute a better approximation to Hks/N, we 
postpone the calculation to a subsequent paper. We mention the new results in Section VI, 
discuss some important issues and indicate some directions for further work. Some of the 
most technical points are presented in the Appendices. 

II. THE SEPARATION OF PHASE SPACE TRAJECTORIES 

Generally, Lyapunov exponents characterize the rate at which infinitesimally close tra- 
jectories in phase space separate or converge in the course of time Q. Here we consider 
a classical system of N identical particles, each of mass m, interacting with short range 
forces in d spatial dimensions. The phase space of such a system is 2dN dimensional and 
any trajectory of the system given by T(t) = (ri(t), Vi(t), r 2 (t), v 2 (t), rW(t), vjv(t)) can be 
specified by giving the Hamiltonian of the system and an initial point, T(0). Here ri(t),Vi(t) 
are the position and velocity of particle i at time t, and we denote a general phase point 
by T. To describe the Lyapunov exponents, we need to consider a bundle of trajectories 
about T(t) which are infinitesimally close to it and which we denote by T(t) + <5r(t). The 
trajectory deviations ST(t) can be expressed in terms of the infinitesimal spatial and velocity 
deviations for each particle as 5T = (Sri(t), 5vi(t), ...,5r N (t),5vN(t)). 

To obtain an expression for the sum of the positive Lyapunov exponents we recall some 
basic ideas of the chaotic dynamics of the phase space trajectories of N particle systems. 
One of the well-known consequences of Liouville's theorem is that the Lebesgue measure of a 
set in the phase space T is constant in time. This implies that if we pick a small set of points 
<5r(0) about an initial phase point T(0), then the measure of this small set fi(5T) will be 
constant in time as the points evolve according to Hamiltonian dynamics. As a consequence 
the sum of all of the Lyapunov exponents must be zero in order that the measure be preserved 
in time. However, if we consider a projection of 5T onto a subspace of lower dimension than 
the full phase space, the measure of this projection in the reduced phase space may change 
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exponentially with time. Suppose then that we project ST(0) onto the velocity subspace. 
Unless this projection is accidentally pathological, its measure in velocity space will grow 
exponentially in time with an exponential given by the sum of all the positive Lyapunov 
exponents of the system (remember that for a time reversal invariant system the Lyapunov 
exponents come in conjugate pairs Aj = —X%, so the larger dN Lyapunov exponents have 
to be > 0). If we denote the projection of ST onto velocity space by SV, and its volume in 
velocity space by ||<5V||, we expect that 

||tfV(t)||~||tfV(0)||exp[t EH (2) 

Ai>0 

asymptotically for large t, where the Lyapunov exponents are denoted by Aj. Before we turn 
to a method for computing the sum of positive Lyapunov exponents, let us write Eq. (|2|) in 
a form that will be useful for our further work. That is, we write the sum of the positive 
Lyapunov exponents formally as a time average. Then by making the assumption that the 
system is ergodic ||, that is, its phase space consists of only one ergodic component, we can 
replace the time average by an ensemble average, as 



EA, = hm - In 
t^oo t 
A,>0 1 



\\SV(0)\\ 



ki^gi (3) 



t^oc t Jo dr \ dr 



where the angular brackets denote an ensemble average, over an appropriate phase space 
ensemble, to be specified further on. If the limit on the left hand side of Eq. (|3D exists and 
is non-zero, the system has positive Lyapunov exponents and is chaotic. Under the above 
assumption about ergodicity the ensemble average will provide an expression for the sum of 
the positive Lyapunov exponents. 

Now we imagine that the dynamics of the system can be characterized as a sequence 
of instantaneous binary collisions connected by free motion of the particles. We neglect 
all effects of collisions where three or more particles are interacting at once, and all effects 
due to the finite duration of the binary collisions. These assumptions are exact for hard 
core potentials but we are neglecting small effects, which are of higher order in a density 
expansion, in the case that the interaction potential is short range but not purely of the 



hard core type. This can be seen by observing that the effects we neglect in such cases are 
of the order of the duration of a collision divided by the mean free time between collisions 
for a typical particle, or of the range of the forces divided by the mean free path length. 
Between collisions the space and velocity variables for all of the particles change with time 

as 

Ti{t) = Vi(t) 

5'nit) = 5v t (t) 
'viit) = 

6iti(t) = 0, (4) 

where we assume, for the time being, that no external forces act on the particles. Since all 
collisions are assumed to be instantaneous binary collisions we can consider the changes in 
these variables at the instant of one such collision between particles k, I say. If we denote the 
variables immediately after the (k, /)-binary collision by primes, then there are no changes 
at the instant of the k, I collision in the quantities fi(t),Vi, Sfi(t), and 6vi(t) for i ^ k,l, but 
for the colliding particles, the changes in these quantities are given by 

= n; = fi 
% = V kl 

if kl = V M ■ v k i 
5R' kl = 5R k r, 5V k \ = 5V k i 
<5r* H = Vki ■ Sfia 
5v kl = Vki ■ 5v k i + ilki ■ Sfia. (5) 

Here R k i = (1/2) (r fc + rj) is the location of the center of mass of particles k, I at the instant 
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of collision, Vki = (1/2) (vk + vi) is the velocity of the center of mass of the two particles, and 
rki = Tk — fi, vm = Vk — vi are the relative position and velocity of particles k, I. Furthermore, 
Vki = 1 — 2<7<r is a reflection operator about a plane orthogonal to the unit vector a, where 



a is in the direction of the line connecting the center of particle I to that of k at the point of 
closest approach of the two particles at their collision and 1 is the unit tensor in d dimensions. 
For a collision to take place we must have (vki ■ cf) < 0. The operator Qki appearing in the 
equation for the change in the relative velocity deviation at collision is given by 
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a a v 



(6) 



[v k i ■ a) 

The above equations for the change in the relative velocity at collision are well known, and 
recently Dellago, Posch, and Hoover have given a derivation for the changes in the spatial 
and velocity deviations at collision, given above, for hard sphere potentials ||. It may look 
surprising at first that Sr 1 ^ ^ 5r k i, because positions do not change in a hard sphere collision. 
The difference is due, however, to the infinitesimal difference in the time of the k, I collision on 
the two nearby trajectories. The results also apply to other short range potentials provided 
one neglects all effects due to the finite duration of a binary collision, including the motion 
of all of the particles during the k, I collision, and provided one identifies a with the unit 
vector along the direction of closest approach between k and /. 

The spatial and velocity deviations are somewhat awkward for calculations and a more 
convenient representation for the rate of separation of trajectories in phase space is obtained 
by supposing that the spatial and velocity deviations are related through a set of matrices 
which we will refer to as radius of curvature, or ROC, matrices. A similar idea was used 
by Sinai in his original formulation of an expression for the Kolmogorov-Sinai entropy of 
a many particle hard sphere system in terms of the ensemble average of the trace of the 
so-called second fundamental operator which depends on the dynamics of the N particles 
in the system 0. Here we define a set of matrices which are more convenient for our 
purposes, but the central idea is certainly due to Sinai. An elementary discussion of the 
second fundamental operator is given by Gaspard and Dorfman ||11|| . We suppose that the 
spatial deviation for particle i at time t, Sfi(t), is related to the velocity deviations of all the 
particles in the system through a set of matrices p i ■ by 
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<^(t) = Pij( f ) ■ MM for i = 1, ...JV. 



(7) 



While this relation couples the spatial deviation of one particle to the velocity deviations 
of all of the particles, we will see below that this coupling is quite manageable, and at low 
density is limited to the coupling of the deviations for a colliding pair of particles. It is 
important to note here that all of the expressions for the sum of the positive Lyapunov 
exponents to be obtained below will depend only on the variables r*j, for i = 1, N, and 
on Pj- for i,j = 1,...,N and not explicitly on the spatial and velocity deviations for the 
particles. 

According to Eq. (Q) we need an expression for the time derivative of ||<5V(t)||. Since 
there are no external forces acting on the system, the volume of the velocity-space projection, 
8V, remains constant whenever the particles are in free motion and only changes with time 
due to binary collisions taking place in the gas. Let us now consider the change in ||5V|| due 
to a binary collision. The velocity deviations after a binary collision of particles k and I are 
related to those before collision by a matrix Am as 



A 



kl 



5v 2 



i 5v N j 



(8) 



where the primed variables denote the velocity deviations immediately after the (k, I) colli- 
sion and the unprimed ones denote the velocity deviations immediately before collision. Of 
course, only the velocity deviations for particles k, I are affected by the collision. Now the 
changes in the velocity deviations for particles k, I are given by Eq. (|5|), which involve the 
spatial deviations for these two particles immediately before collision. These, in turn, are 
related to the velocity deviations of all the particles through the ROC matrices defined in 
Eq. (^). Putting these equations together and noting that the volumes in velocity space 
before and after the (k,l) collision are related by 



\8V\\ = | det Afci| ||5V||, 



(9) 



we find that we will need the determinant of A k i which is given by 



det Am = det 



U fci + -fl kl ■ {p kk + p u - p kl - p lk ) 



(10) 



Besides Eqs. (§) and (|7|) we have used here the relation, 5V kl = 6V k i as well as the fact that 
the velocity deviations of all the other particles are unaffected by the (k, I) collision. In Eq. 
(0) above, we have expressed the sum of the positive Lyapunov exponents in terms of the 
ensemble average of the time derivative of In ||5V||. Since this quantity only changes at the 
instant of a binary collision, we need to (a) look at the binary collisions taking place in the 
gas at some instant; (b) calculate the change in In \\SV\\ at each of them; and (c) take the 
average over an appropriate stationary ensemble distribution. Since we have a total of N 
particles in the system we need to count the possible contributions from each of N(N — l)/2 
possible binary collisions in the gas. Thus we find that the sum of the positive Lyapunov 
exponents is given by 

12 X i =12 I da\v kr a\e(-v kr a)S(f kl - aa)ln\ det A k i\). (11) 

Ai>0 k<l ^ J ' 

Here the step function 0(x) = 1 for x > and zero otherwise. We have integrated over all 
possible values of a for each collision, included the fact that the centers of the particles must 
be separated by a distance of a at collision, and included that rate at which collisions take 
place in the factor \v k i • a\ in order to get the proper expression for the time derivative of 
the velocity space volume. This expression is exact for hard-sphere systems in any number 
of dimensions, and can be used, with appropriate modifications to the differential scattering 
cross section, as a good approximation for other short range potentials. 

To proceed further we must explain how the ensemble average is to be computed. We 
suppose that it is possible to find a stationary distribution function which describes not 
only the space and velocity variables, fi, Vi, for each particle, but also the distribution of the 
N 2 ROC matrices. While this seems especially complicated, it will turn out that only some 
reduced one and two particle distribution functions are actually needed for our computations. 
We assume the existence of a distribution function J r N (xi, xjv, Pn, Pvii ••■> Pnn)i which is 
normalized as 
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J dx 1 ...dx N dp 11 ...dp NN J 7 N = 1. (12) 

here Xi = (r*j, Vi) and the integrations are over all components of position and velocity for each 
particle as well as over all the elements of each of the ROC matrices. A simplification occurs 
immediately since the ensemble average is over functions of two particles. If we assume that 
Tn is symmetric in its variables we can express the sum of the positive Lyapunov exponents 
in terms of a two particle distribution function T<i{x\, x 2 , Pn, Pm P211 P22X as 

A,; = / dxidx2dp n ...dp 22 / da\vi 2 ■ <r|©(— Vi 2 ■ &)^2^ifi2 — hi I det A 12 |. (13) 

\,>o 2 J J 

Here the two particle distribution function is defined by 

N\ r 

J r 2 (xi,x 2 , p u , Pi2> P2i ? P22) = , N _ 2 y J dx 3 ...dx N H dp^Txixx, p NN ), (14) 

where the prime on the product means that the elements of the matrices p n , p 12 , p 2 i, P22 are 
not to be included among the integration variables. This is the main result of this section 
- an expression for the sums of the positive Lyapunov exponents in terms of an extended, 
stationary two particle distribution function. In the next section we shall derive the BBGKY 
hierarchy equations that are needed to determine this function. 



III. THE BBGKY HIERARCHY EQUATIONS 

For hard sphere systems one can describe the time evolution of any phase space quantity 
in terms of a time evolution or "pseudo-Liouville" operator, C + (N), such that any dynamical 
quantity A(T(t)) that does not have an explicit time dependence changes with time as [|12| -[T4]| 

A(T(t))=exp[tC + }A(T(0)). (15) 

Here the subscript + indicates that we are considering the forward time evolution, t > 
0, of the system, since the operators differ slightly for evolution with t < 0, due to the 
different ways they are defined for nonphysical initial configurations. The phase point T(t) 
is physically meaningful only if no pair of particles is separated by a distance less than a, 
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but even for unphysical phase points the time evolution operator is well-defined. To be sure, 
for obtaining physically meaningful results one has to multiply expressions of the form of the 
right hand side of Eq, ( |T5| ) by a weighting function containing the factor (r) defined to be 
zero for overlapping configurations and unity for non-overlapping configurations. Here we 
extend the definition of the pseudo-Liouville operator to include the time dependence of the 
ROC matrices in addition to the positions and velocities of all of the particles. We will then 
discuss the adjoint of this operator which is used to obtain the BBGKY equations for the 
extended distribution functions T n for n = 1,2, .... Since the method for constructing the 
pseudo-Liouville operator has been discussed in the literature at some length, and since the 
calculation here is an obvious extension of this method, we merely outline the calculation 
here, following the procedure in Ref. ||14|| . 

Suppose we had a system with only one particle. Then the time evolution of some 
function A(fi, vi, p n ) is trivial: 

A(fi, vi, p n ,t) = A{ri + v x t, v u p n + It) = exp[t£ (l)]A(rL, v x , p n ) (16) 

where 

d d d 

u 'l a =l °Pll,aa 

Here we assume that there are no external forces acting on the system, the derivatives 
are taken with respect to the spatial coordinates and with respect to the diagonal ele- 
ments of the ROC matrix. This latter condition follows from the observation that un- 
der free particle motion, only the diagonal parts of the ROC matrices change with time, 
growing linearly as indicated above. Now consider the less trivial case of two particles. 
Let A(xi, x 2 , Pu, p u , P211 P22) be some function of the indicated variables and consider the 
quantity I t A given by 



XtA 



S t (l, 2) - <S°(1, 2)1 A(x u x 2 , p n , p 22 ) (18) 



where «S t (l, 2) replaces all of the dynamical variables by their values at a time t later, 
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obtained by following the motion of particles 1 and 2 forward in time, and Sf(l, 2) is a free 
streaming operator that acts on A to produce 

2)A(ri, tfi, r 2 , r7 2 , Pn, Pi 2 , P 2 i, P22) 
= A(fi + wit, vi, f 2 + w 2 t, tf 2 , p n + It, p 12 , p 21 , p 22 + It) 

= exp[t£ (l,2)]A(fi,...,p 22 ). (19) 

where £o(l,2) = £o(l) + £o(2). By considering the two particle dynamics and by using 
the method of Ref. ||14]| , we can express in terms of the free streaming operators and a 



binary collision operator 7+ 2 ^(l,2) as 

l t A = f drS° T (l, 2)Tj 2) (l, 2)S°_ t (1, 2)A(f u p 22 ), (20) 

where the the binary collision operator is given by 

7^(1, 2) = a d - x J da\v 12 ■ &\G(-v 12 ■ &)5(f 12 - aa)(V a (l, 2) - 1). (21) 

In Eq. ( PH ) V a (l, 2) is a substitution operator which replaces Vi, v 2 , p n , p 12 , p 21 , p 22 by their 
values immediately after a (1, 2) collision. These values are given by Eq. (^|) as 

A = #l - (#12 • 0")<7 

"4 = ^2 + (#12 ■ o-)o" (22) 

for the velocities of the two colliding particles, and by the following equations for the ROC 
matrices, derived in Appendix A, 

Pa + Pb ■ Ol2 • P d = P a , 

p' b ■ U12 + p' b ■ flu ■ P c = P&, 

Pc • Ui 2 + p^ ■ Oxa ■ p c = U12 ■ p c , 

Pd + Pc' O12 ■ P d = U12 ■ p d , (23) 
where the new ROC matrices are linear combinations of the given by 
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Pa = g(Pll + Pl2 + P21 + P22) 
Pb = \{Pll ~ Pl2 + P2I - P22) 
Pe = 2VP1I ~ Pl2 - P21 + P22) 

Pd = (P11 + P12- P21- P 22 )- ( 24 ) 

As usual, the primes denote values of the quantities immediately after a collision. 

The superscript (2) on the binary collision operator denotes the condition that we are 
only keeping two-particle contributions to the ROC matrix expansions for 5ri and 5f 2 . In 
principle we will have to keep many particle contributions to the ROC matrix expansions 
of the spatial deviations which will force us to define a set of binary collision operators, 
Tf n) (l,2|3,4,...,n), for n = 2, ...,7V, even for a collision which only involves particles 1 and 
2, since the ROC matrices must be "updated" at each collision. While this causes some 
complications in the formalism, these complications can be dealt with along the lines used 
in the statistical mechanics of systems with more than two body forces. Appendix B outlines 
the cluster expansion technique needed to handle these extended binary collision operators, 
but for most of the purposes of this paper, this cluster expansion will not be needed. We 
mention here that unless there was some previous dynamical connection of particle k, say, 
with either particle 1 or 2, the ROC matrices p lk , p 2k , p kl , p k2 will all be zero before and 
after the (1,2) collision. 

Now that we have defined a binary collision operator that acts on functions of dynam- 
ical variables, we take its adjoint to get an operator that acts on distribution functions. 
That is, consider an average value of S t (l, 2)A(xi, p 22 ) taken with respect to some initial 
distribution function jF 2 (a;i, p 22 , t = 0) which we assume contains a factor W(l,2), the 
two particle overlap function. The adjoint of the two-particle time displacement operator is 
defined by 

J dx 1 dx 2 dp 11 ...dp 22 J 7 2 (xi, p 22 , t = 0)«S t (l, 2)A(x 1: p 22 ) 

= J dx 1 dx 2 dp u ...dp 22 A(x 1 , x 2 , p n , p 22 )«S_ t (l, 2)J r 2 (x 1 , p 22 , t = 0), (25) 
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where the "bar" denotes an adjoint, and the subscript notation on the adjoint of the stream- 
ing operator calls attention to the fact that while dynamical variables "move forward" in 
time, distribution functions " move backward" in time. 

After some elementary calculations one finds that there is a binary collision operator 
representation of the adjoint streaming operator of the form 

2) = S°_ t (l, 2) + drS°_ T (l, 2)fi 2) (l, 2)S°_ (t _ T) , (26) 

where 

S£ t = exp[-t(£o(l,2))] and (27) 
Ti 2) (l,2) = o d ~ 1 J da\v 12 ■ a\Q(v 12 ■ a)x 

5(^2 -era) f dp' u dp' 12 dp' 21 dp' 22 J] S(p lJ -p lJ (p' n ,...,p' 22 ))V' a (l,2)~S(f 12 + aa) .(28) 

Here V' a (l,2) is a substitution operator that replaces ROC matrices to its right by the 
corresponding primed matrices, and velocities by the restituting values, namely those which 
for a given a produce vi,v 2 after a binary collision. The delta functions in the ROC variables 
require that the primed ROC matrices be the restituting ones, i.e., those which produce the 
unprimed ROC matrices after a binary collision. These restituting ROC matrices are found 
by interchanging the primed and non-primed matrices in Eq. (0), and solving for the 
primed matrices in terms of the unprimed ones. 
The operator 7+ (1,2) has the property that 

Tj 2 \l,2)S T Tj 2 \l,2) = (29) 

since it is impossible for particles 1 and 2 to collide more than once without the intervention 

— ( 2 \ 

of a third particle. A similar identity holds for the 7_ (1, 2) operator, so that we may write 
the time displacement operators <St(l, 2) and «S_t(l, 2) in a pseudo Liouville operator form 

<S t (l, 2) = exp t[£ (l, 2) + T?\l, 2)], (30) 

whenever the time displacement operator acts on dynamical variables and a function W(l, 2) 
appears to the left of the time displacement operator, and 

15 



S- t (l, 2) = exp -t[C (l, 2) - Ti 2) (l, 2)], (31) 

whenever the time displacement operator acts on distribution functions which vanish for 
r 12 < a. 

This analysis of the two particle time displacement operators can be readily extended to 
the N particle case, and it is this extension which allows us to derive a set of BBGKY hier- 
archy equations for the s particle distribution functions defined by a natural generalization 



ofEq. PP 



T s {x\, x„ p n , p ls , p 21 , p ss , t) = ; J dx s + 1 ...dx N 'Y^dp i:j J : N {x 1: p NN , t) 



(32) 



where we have included a possible time dependence of the distribution functions. The 
normalization used here is such that the successive distribution functions, T s are proportional 
to (n<T d ) s , that is to the s-th power of the reduced density of the gas. 

To derive the BBGKY equations we start with the N- particle pseudo-Liouville operator 
and write the extended Liouville equation as 

d 



^-+jC (l,2,...,N)-J2^ N \iJ\h-,N) 



Tn = 0, (33) 



where 2, N) = J2i=i^-o{i) and the N-particle binary collision operators are defined 
so as to include the effects of the collision of particles i and j on all of the ROC matrices. 
We hasten to add that unless there was some previous dynamical connection of particle k, 
say, with either particle i or j, the ROC matrices p ik , pj k , p k j, p ki will all be zero before and 
after the collision. This greatly simplifies the use of these complicated looking binary 
collision operators. In Appendix B we show that the iV-particle binary collision operators 
can be expressed in terms of a cluster expansion of the form 

N 

T± N) {1, 2|3, ...,N) = Ti 2) (l, 2) + 2|3, n) (34) 

n=3 

where the n-particle cluster functions are defined recursively in Appendix B. Here we give 
only the first one, 
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Wi 3) (l, 2|3) = ri 3) (l, 2|3) - Ti 2) (l, 2). (35) 

If we now integrate over the pseudo-Liouville equation, Eq. ([53|), and use Eqs. Q3"2"D, and 
([35]), we obtain the BBGKY hierarchy equations. In the following we will need only the first 
one 



^l(^l.Pll) 



N 



J dx 2 dp 12 dp 2l dp 22 Tl 2 \l,2)J r 2(x 1 iX2,Pn,-"p22) + l>2 J dx ? ,...dp nn U { ? ) {I, 2|3, n) T n (36) 

n=3 

We have not given the explicit form of the remaining hierarchy equations, but they can easily 
be obtained, if so desired. When these equations are integrated over all of the elements of the 
ROC matrices, then the usual BBGKY equations are recovered. The extensions given here 
can be used for a systematic determination of sums of Lyapunov exponents for dilute and 
moderately dense gases. In addition they allow for non-systematic approximations, similar 
to Enskog's equation for hard sphere systems |10|JT7|] , which may possibly be used even for 
systems in glassy states. 

In the next section we truncate the hierarchy at the first equation to obtain an extended 
Boltzmann equation for the sum of the positive Lyapunov exponents at low density. 



IV. THE EXTENDED BOLTZMANN EQUATION 

In the previous section we outlined the derivation of the full set of the extended BBGKY 
hierarchy equations. Their solution is a formidable task, which has not been fulfilled even 
for the usual hierarchy equations. Nevertheless substantial progress has been made toward 
approximately solving the usual equations, leading to a more basic understanding and im- 
portant modification of the Enskog theory of dense hard sphere systems, to explanations 
of long time tail phenomena in dense fluids, and to theories of glasses, among other results 
T0| , [r5|| . Here we will not attempt to carry out a similar analysis for the extended equations 



- leaving that for another place - but will restrict our attention to the theory for very low 
densities. 
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We note from Eq. ([13]) that the 2-particle distribution function JF 2 , is needed in order 
to determine the sum of the Lyapunov exponents. Specifically, we need this function at 
the instant of a binary collision of two particles (1, 2), say. From experience with the usual 
hierarchy equations we propose that a useful expression for T 2 at the instant just before 
collision in Eq. ( |T3"D is given by 

•T^Oi, x 2 , p n , p 12 , p 2 i, Pmt) = J r i(xi,p 11 ,t)J'i(x2,p 2 2,t)5(p 12 )5(p 21 ) 

+G 2 (x 1 ,x 2 ,p U} ...,p 22j t), (37) 

where we have introduced a new cluster function Q 2 . In the first term on the right hand 
side of Eq. ( |3"7| ) we assume that particles 1 and 2 are totally uncorrelated at the instant 
of the collision. This not only leads to the factorization of the pair distribution function, 



but also to p 12 = p 21 = 0, giving rise to the delta functions in Eq. (p7|) . Any effects 
of a possible correlation between the particles are put in Q 2 . The Boltzmann equation 
approximation, which we will use here, is to neglect the contribution of the cluster function 
Q 2 in the expression for the sum of the Lyapunov exponents. We then need to determine 
the one-particle distribution function T\, which satisfies the first hierarchy equation, Eq. 
(j36|), where the right hand side depends upon two and higher particle distribution functions. 
Again, we use the Boltzmann approximation for JF 2 , and neglect the contributions from the 
higher order cluster operators and distribution functions on the right hand side of Eq. fl36|). 
These approximations can only be justified after one has carefully examined the terms that 
are neglected, and shown that they are of higher order in the density. We leave this for a 
further publication, and explore here the effects of these approximations. 

Under the above Boltzmann approximations we have an expression for the sum of the 
positive Lyapunov exponents for low densities following from Eqs. @, @, (|T0|) as 

— ~7j— J dx 1 dx 2 dp u dp 22 J da\v 12 ■ a\Q(— v 12 ■ a)5(f 12 - a a) x 

Pn) ^1(^2, P 22 )- ( 38 ) 



U12 + -Oia • (Pn + p 22 ) 



In I det 

Here the one particle distribution functions obey, in this approximation, a steady-state, 



extended Boltzmann equation of the form for i = 1,2, 

This equation has to be supplemented by a boundary condition at p u = 0. For pa < 0, T\ 
has to vanish since semi-convex scatterers never produce negative radii of curvature if those 
were not present initially. The boundary conditions at infinity follow from the requirement 
that the distribution functions, when integrated over the ROC variables, reduce to the usual 
space and velocity distribution functions. We use a steady state distribution function since 
we are interested in the sum of the Lyapunov exponents in either equilibrium states, or non- 
equilibrium steady states. Here we consider an equilibrium state, and look for solutions of Eq. 
(|39|) which, upon integration over the elements of p u reduce to the usual Maxwell-Boltzmann 
equilibrium distribution function. Once we have determined T\{xi, p u ) for i = 1, 2 by solving 
Eq. (|39"D , we can compute the sum of the positive Lyapunov exponents using Eq. (|38|). 

V. THE SIMPLE COLLISION DAMPING APPROXIMATION 

In our previous work on the Lyapunov spectrum for the random Lorentz gas at low 
densities, we obtained an extended Lorentz-Boltzmann equation [^]. This equation could be 
solved under equilibrium and non-equilibrium steady state conditions to provide expressions 
for the Lyapunov spectrum for the Lorentz gas which agree very well with the computer 
simulation studies of Posch and Dellago. There we found that the dominant contribution to 
the Lyapunov spectrum comes from the collisional damping of the appropriate distribution 
functions for the elements of the ROC matrices, obtained either by solving the extended 
Lorentz-Boltzmann equation for the distribution of the ROC matrix elements, or from a 
more heuristic kinetic theory analysis which shows that the distribution of the ROC matrix 
elements is closely related to the free path distribution for particles in the gas. A similar 
free path analysis was applied to the sum of positive Lyapunov exponents for a dilute gas in 
equilibrium. There it was found that the analysis provides an expression for this sum which 
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is of the form given by Eq. (p]) with excellent agreement with simulations for the coefficient 
a, but less than satisfactory agreement for b in Eq. ([I]). Here we show that the free path 
results can be easily recovered from Eqs. (^) and fl39|), by making a simple approximation 
in the expression for the binary collision operators appearing in Eq. (|39|). A better analysis 
of this equation leads to improved values for b, while leaving the values of a unchanged. 
This improved analysis is somewhat lengthy and will presented elsewhere. Here, we show 
how the previous results are recovered. 

We begin by noticing that the expression, Eq. (28), for the binary collision operator 
T_ (i, j) consists of two parts, a restituting term which replaces variables by their restituting 
values, denoted by primed variables, and a direct part which accounts for the loss of particles 
with the desired values of the variables due to collisions. As the above remarks suggest, a 
simple approximation to the extended Boltzmann equation, which incorporates the effect 
of collisional damping of the distribution of the ROC matrix elements, may be obtained by 
assuming the restituting part of the binary collision operator on the right hand side of Eq. 
is proportional to to <5(p 11 ). Such an approximation leads immediately to 

£ Q .Fi (xi, p u ) « -a^ 1 J dx 3 dp 33 da\v i3 ■ ■ a)5(f i3 + oa)Fi(xi, p H ) Fi{x 3 , p 33 ). (40) 

Now we use the fact that the system is assumed to be in equilibrium. Then we can assume 
that T\{xi, pa) does not depend upon r«, and that when one integrates over the ROC vari- 
ables in the one particle distribution function, one obtains the Maxwell-Boltzmann velocity 
distribution function. That is 

J dp 33 F x (v 3 , p 33 ) = mp(v 3 ), (41) 
where n = N/V is the number density of the gas, and 




where f3 = (ksT) , with T the thermodynamic temperature of the gas, and ks is Boltz- 
mann's constant. Thus we are led to a simple, approximate, differential equation for 
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d q 

Yl o Fi( x i> Pa) ~ Pa), (43) 

a=l @Pii,aa 

where is the equilibrium, low density collision frequency for a particle with velocity v% 
given by 

u(vi) = na d - 1 J dv 3 J da\v 13 ■ a\G(-v 13 ■ a)(p{v 3 ). (44) 

In order to complete the calculation of the sum of the Lyapunov exponents, we return 
to Eq. (H|) and examine the determinant in more detail. Since the operator U12 is its own 
inverse and the magnitude of its determinant is unity, it follows that 

|det[U 12 + 12 -(p n + p 22 )/2]| = |det[l + r 12 -(p n + p 22 )/2]| where (45) 



a 



v 



Vi 2 a + avu - icr ■ Vi 2 )l - 77 — =r-r<7<7 



(46) 



[a ■ v 12J 

It is easily seen that the operator T 12 is proportional to a projection operator onto a subspace 
orthogonal to the relative velocity vector {T 12 . We consider the two and three dimensional 
cases separately since a bit of geometry is required now to complete the calculation. 

We begin with the case of two dimensions. Then the determinant we need is found to be 



2 


V12 




a 


COS0 



det[l + r 12 ■ (p n + p 22 )/2] = 1 + -p! l -p ±± where (47) 



v 12 ■ a = \v 12 \ cos0; and p ±± = ^(um • (pn + P22) • "izl- (48) 

Here the subscript lona two dimensional unit vector vu,± denotes a unit vector orthogonal 
Vu- From Eq. ( |4lj| ) we see that we need the distribution for only one of the diagonal com- 
ponents of the ROC matrices p n and p 22 in a representation based on the coordinate frame 
defined by , 0i2,'0i2±- From Eq. ([4"3"D, integrated over the other component, this distribution 
function, properly normalized, is easily found to be 

Fifa, pu,_L±) = nv(vi)ip(vi) expf-z/^p^xj. (49) 

When one uses this expression in Eq. and calculates the quantity J2\i>o K/N,the 

KS entropy per particle for a dilute gas of hard disks, one recovers exactly the expression 
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obtained in Ref. 0, given in Eq. (p]) above. We note that the first term on the right hand 
side of Eq. (0), can be neglected here since the next term is on the order of the free 
time between collisions divided by the time needed to travel a molecular diameter which is 
of order (r/\vi2\- For a dilute gas this ratio is much larger than 1 and the error made in 
neglecting the first term is higher order in the density. 

The three dimensional case is treated in a very similar way. In this case the space or- 
thogonal to V12 is two dimensional, and the corresponding operator Tn projects the matrices 
Pn, p 22 onto this subspace. Thus, instead of needing the nine components of each of these 
matrices, we need only four. We then find, for three dimensions 

2 \&& - Pf±Pf±] (50) 



det 
where 



1 + • -(P11 + P22) 



2k 2 | 



a 



Jj _ . (PU + P22) _ Pll±± + P22±± 



P±± — V 12,± 2 V 12,± — 2" ' ^ ' 

and the unit vectors t>i2, v^ ± , v^x f° rm an orthogonal coordinate basis in three dimensions. 
Of the four matrix elements needed for each of the ROC matrices only the two diagonal 
components (in superscript) are important here. This follows from the observation that, 
according to Eq. (|17D , between collisions the diagonal components grow linearly with time, 
while the off-diagonal elements remain constant, changing only at collisions. Since the 
time between collisions scales inversely with the density, only the diagonal components are 
important at low density for obtaining the dominant contribution to the sum of the Lyapunov 
exponents. We then find that 



a 2 r 

hxg/N = — J dvidv2dp\\ ±±dp^ tl _j_da\vi2 • o"|0(— ^12 ■ o) In 



1^12 I (Pn ,_LL +P22,±± 



a 



The single particle distribution functions appearing in Eq. ( |52|) have the same form as those 
in Eq. ( pEUp with the three dimensional form of the collision frequency used in place of the 
two dimensional form. This too is the same result as that obtained in Ref. || 
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VI. CONCLUSION 



We have outlined here a method that allows a systematic expansion of the sum of the 
positive Lyapunov exponents for function of the gas density. For a gas in equilib- 

rium, this leads to an expression for the KS entropy which, until recently ||, had been the 
subject of much speculation but with no explicit analytic results ||16|| . We have formulated 
the theory for a gas of hard spheres but it can be extended to gases which interact with 
other short range potentials by extending the definition of the binary collision operators to 
such potentials. Such an extension will almost certainly ignore effects that take place on the 
time scale of the duration of a binary collision as well as effects of genuine many-particle 
collisions. Nevertheless, such an extension would provide useful, if not completely accurate 
results. A similar situation obtains in the kinetic theory for transport coefficients and time 
correlation functions where hard sphere results can be extended to other potentials with 



similar approximations [10]. A simple way to extend the low density results to higher den- 
sities, for hard sphere systems and other systems with short range potentials, is to use the 
Enskog approximation when computing the two-particle distribution function |fT0|JT7| ; p~8[1 . In 
this approximation one replaces the two particle function at the instant of a collision by 
the product of single particle functions as done here, but multiplies this product by the 
pair distribution function for two hard spheres at contact, the so-called Enskog ^-factor. 
This approximation leads to useful approximations for the transport coefficients and time 
correlation functions of dense gases 

There are now several clear directions for future work: 

1) The low density results as derived here must be improved by avoiding the approxima- 
tion that the restituting collisions are proportional to <5(pjj), as indicated in Section IV above. 
Some progress has been made in this direction. One of us (JRD) found, in an approximate 
calculation of the restituting contributions, that the coefficient a in Eq. ([!]) is unchanged 
as expected, but there is an addition to b which is about In 4 in two dimensions, changing 
b from 0.2 to ~ 1.6, compared to the computer result of b ~ 1.35, In three dimensions, an 
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additional contribution to b of approximately In 3 is found, changing the previous result from 



b = 0.56 to b ~ 1.66, compared to the computer result of b ~ 1.34 [21]]. These results are 
now much closer to the computer results and an encouraging sign of the correctness of this 
approach. It is worth remarking here that the important role of the restituting collisions is 
due to the "semi-dispersing" nature of hard sphere systems, which is quite different from 
the Lorentz gas system which is a purely dispersing system. In a collision between particles 
k and I not all of the components of the ROC matrices p kk and p u will be reduced to values 
on the order of 0"(m/fcsT) 1,/2 ), as is the case for the Lorentz gas. This will be explained in 
a further publication |f22fl . 

2) These methods and results must be extended to other situations where sums of Lya- 
punov exponents play an important role. Such situations include trajectories on the fractal 



repeller in an open system, as discussed by Dorfman and Gaspard [23]. In such systems the 



sum of the positive Lyapunov exponents together with the KS entropy (which does not sat- 
isfy Pesin's theorem in this case) determines the transport coefficients of the system. Other 
situations of interest occur in systems subject to external driving forces in combination with 
Gaussian thermostats. In such systems the Liouville measure is not preserved. There is an 
interesting connection between the sums of all the Lyapunov exponents and transport phe- 
nomena such as entropy production and transport coefficients for such systems, which has 
been discussed in detail in the literature [f24]| . The methods discussed here can be extended 



to treat such thermostatted systems, also, but one has to calculate both the sums of all the 
positive and all the negative Lyapunov exponents. 

3) The method used here can be extended to higher densities by including the suc- 
cessively higher hierarchy equations and correlation functions. In addition one can make 
non-systematic high density approximations in the spirit of the Enskog equation. Among 
other reasons it is important to look at higher density effects since one can then make some 
connections between dynamical systems theory and long time memory effects in a fluid. 
This will be important not only for a deeper understanding of non-equilibrium phenomena 
in general, but also for answering some of the questions that arise in the theory of glasses, 
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particularly those associated with the ergodic properties of such systems . 
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APPENDIX A: DERIVATION OF EQ. (j24|) 



We begin by noting that for a system of two particles, the spatial deviation vectors are 
expressed in terms of ROC matrices as 

Sri = Pn ■ Sv u + p 12 • 5v 2 , 

Sr 2 = p 2 i • Sv n + p 22 ■ 5v 2 . (Al) 

From this it follows that 

5R U = p a ■ 5V 12 + p b ■ 5v 12 , 

Sr n = p c ■ Sv 12 + p d ■ 5V 12 , (A2) 



where p a ,...,p d are defined in Eq. (|24]). Now consider the changes in 5R\ 2 and 5r i2 at 



collision as given by Eq. (5). The center of mass equation becomes - the primed variable 
denote values after collision - 

Pa ■ SV[ 2 + p' b Sv 12 = p a ■ 5V 12 + p b ■ 5v 12 . (A3) 

We now insert the expressions SV[ 2 = 5V\ 2 and and the appropriate one for 5v[ 2 given by 

5v 12 = U12 • V12 + O12 • [p c ■ Sv 12 + p d ■ 5V 12 \, (A4) 

into Eq. (|A3|), to obtain 
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Pa ■ $Vi2 + p' b ■ [U12 • Sv 12 + ^12 • (Pc * ^13 + Pd * ^12)! = Pa * <^12 + Pfe * &?12- (A5) 



We can now equate coefficients of SV 12 and of Sv^ on each side of Eq. ( |A3f) since the relative 
and center of mass velocities are independent of each other, and there are no orthogonality 
or normalization contraints on the components of each of these velocity deviation vectors 
which would make some of the components of each of these vectors dependent on the others. 
One therefore finds the first two lines of Eq. (|23|). A similar treatment of the collision 
equation for <5r 12 leads to the third and fourth lines of Eq. fl23|). 

APPENDIX B: CLUSTER EXPANSIONS OF THE BINARY COLLISION 

OPERATORS 

Cluster expansion methods have proven to be of great value in both the equilibrium 
and the non-equilibrium statistical mechanics of fluid systems |T0| , p6| - p8[| . They allow us to 



express functions defined for a system of N particles in terms of functions defined for systems 
of fewer numbers of particles. In equilibrium fluids, for example, one can obtain the virial 
expansions of thermodynamic functions using Ursell cluster expansions of the iV-particle 
Boltzmann factor exp[— (3H^\. Here we apply the Ursell method to the N particle binary 
collision operator 7j (1, 2|3, 4, N). The idea is to choose the first term in the expansion 
as a two particle function, with the next and higher order terms successively correcting for 
the errors made in the earlier terms. For example, for three particles, we would write 

Tj 3) (l,2|3) =Ti 2) (l,2)+ZYi 3) (l,2|3), (Bl) 

where, obviously, 

Wi 3) (l,2|3) = Ti 3) (l,2|3) -Ti 2) (l,2). (B2) 

This procedure can be repeated for four particles to define a four particle cluster operator 
W_ (1, 2 1 3, 4) through the equation 

Ti 4) (l,2|3,4) = Ti 2) (l,2)+wi 3) (l,2|3)+wi 3) (l,2|4)+wi 4) (l,2|3,4). (B3) 
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Extending this idea to N particles, we find 

TV 

Tl N) {l, 2|3, 4, N) = T (2) (l, 2) + £ Wi 3) (l, 2 1 Jfe) + ]T U^\l, 2| fc, /) + ••• . (B4) 

fe=3 k<l 

This cluster expansion of the binary collision operator is used to obtain Eq. fl36|). Similar 
cluster expansions, with simple modifications, can be used for the N particle distribution 
functions and pseudo Liouville operators. 
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